#include <arm_neon.h>
#include <sys/time.h>
#include <iostream>
#include <random>

using namespace std;

#define N 4096

float A[N][N];

void createMatrix() {
    srand(20000214);
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            A[i][j] = rand() % 100;
            // cout << A[i][j] << " ";
        }
        // cout << endl;
    }
    // cout << endl;
}

void LU() {
    for (int k = 0; k < N; k++) {
        for (int j = k + 1; j < N; j++) {
            A[k][j] /= A[k][k];
        }
        A[k][k] = 1;
        for (int i = k + 1; i < N; i++) {
            for (int j = k + 1; j < N; j++) {
                A[i][j] -= A[i][k] * A[k][j];
            }
            A[i][k] = 0;
        }
    }
}

void Neon_LU() {
    float32x4_t t1, t2, t3, t4;
    for (int k = 0; k < N; k++) {
        float kk = A[k][k];
        t1 = vdupq_n_f32(kk);
        for (int j = N - 4; j >= k; j -= 4) {
            t2 = vld1q_f32(A[k] + j);
            t3 = vdivq_f32(t2, t1);
            vst1q_f32(A[k] + j, t3);
        }
        for (int j = k; j % 4 != N % 4; j++) {
            A[k][j] /= kk;
        }
        for (int i = k + 1; i < N; i++) {
            float ik = A[i][k];
            t1 = vdupq_n_f32(ik);
            for (int j = N - 4; j > k; j -= 4) {
                t2 = vld1q_f32(A[i] + j);
                t3 = vld1q_f32(A[k] + j);
                t4 = vsubq_f32(t2, vmulq_f32(t1, t3));
                vst1q_f32(A[i] + j, t4);
            }
            for (int j = k + 1; j % 4 != N % 4; j++) {
                A[i][j] -= A[i][k] * A[k][j];
            }
            A[i][k] = 0;
        }
    }
}

int main() {
    struct timeval t1, t2;
    createMatrix();
    gettimeofday(&t1, NULL);
    LU();
    gettimeofday(&t2, NULL);
    double deltaT = (t2.tv_sec - t1.tv_sec) * 1000000 + t2.tv_usec - t1.tv_usec;
    cout << deltaT / 1000 << endl;
    /*for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++)
            cout << A[i][j] << " ";
        cout << endl;
    }*/
    gettimeofday(&t1, NULL);
    Neon_LU();
    gettimeofday(&t2, NULL);
    deltaT = (t2.tv_sec - t1.tv_sec) * 1000000 + t2.tv_usec - t1.tv_usec;
    cout << deltaT / 1000 << endl;
    /*for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++)
            cout << A[i][j] << " ";
        cout << endl;
    }*/
}
